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This paper examines several of the approximations or models Involved In 
the development of a numerical combustor flow code. In the first section, the 
Importance of numerical accuracy Is Illustrated, and the Impact that Improved- 
accuracy schemes have on slowing convergence Is demonstrated. Solution algo- 
rithms that can speed convergence are discussed and some performance features 
of these algorithms are Illustrated. A sample calculation displaying the 
Importance of boundary conditions on a three-dimensional numerical prediction 
Is presented. The Inaccuracy of a current turbulence model In highly turbu- 
lent (nonequilibrium) regions Is described. Finally, the surprisingly good 
performance of a six-flux model In describing radiation heat transfer Is dis- 
played. In all the areas examined, continued research Is still needed, but 
valuable engineering tools are available today. 

INTRODUCTION 

Three-dimensional combustor calculations Involve detailed modeling of 
several Important physical processes. Airflow, chemical reactions, fuel 
sprays, and turbulence are just a few of the physical processes that must be 
described. Many of these processes occur on both a molecular and a macro- 
scopic scale. To exactly describe these processes numerically, one must 
resolve these scales on a computational mesh. And this Is clearly beyond 
current computational resources. To make the computational task tractable, 
we Introduced modeling assumptions. These modeling assumptions limit the 
generality of the computational flow code, but It Is hoped that the dominant 
physics remain correctly represented. 

Modeling assumptions are only the first limit of generality Introduced 
when developing a combustor flow code. A further limit Is Introduced by the 
need to approximate the modeled equations before they are solved numerically. 
This approximation process can significantly affect the accuracy of a model 
prediction. 

In view of all these factors, a combustor designer cannot be expected to 
fully embrace a computer model prediction. If a clear distinction between 
modeling errors and errors Introduced by the numerical solution algorithm can- 
not be made, then the designer Is left with a very unreliable computational 
tool . 


This paper examines some of the compromises made when developing combus- 
tor flow codes and how these compromises affect the accuracy of turbulent flow 
calculations. To assess the balance of the modeling compromises, we compare a 
number of sample calculations with experimental data. This paper examines the 
five areas of numerical accuracy, solution algorithms, Inlet boundary condi- 
tions, turbulence models, and radiative heat transfer. 
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The NUMERICAL ACCURACY section focuses on the development of more accu- 
rate numerical methods to be used In combustor flow codes. Upwind differenc- 
ing, which Is currently used in these flow codes, Introduces an appreciable 
error (or numerical diffusion) Into the calculation. This error may be of 
such a large magnitude that It obscures the turbulence model used In the cal- 
culation. A series of calculations are Illustrated which demonstrate the 
accuracy of a variety of differencing schemes. 

An Important aspect of Improved accuracy Is the effect that the differ- 
encing scheme has on the rate of convergence. The Improved accuracy schemes 
all appear to require more CPU time to converge. For this reason, the next 
section discusses solution algorithms. SIMPLE ( seml-lmpl 1 c 1 t pressure-linked 
equations) Is one of the most widely used solution algorithms for solving the 
steady-state form of the Navler-Stokes equations. Although this scheme has 
proven to be quite effective. Its convergence rate can be Improved. This 
paper focuses on two approaches which accelerate convergence by performing 
corrections that Improve the Iterative agreement with continuity. These 
alternate schemes are Illustrated In a series of calculations. 

A third section examines the Importance of Inlet boundary conditions. An 
Illustrative example Is displayed. 

The fourth section discusses turbulence models and reaction closures. 

Here the more pragmatic approaches to calculating turbulent reactions are 
Illustrated. An eddy-breakup model and a PDF (probability density function) 
method are described and compared. 

Finally, the fifth section examines the accuracy of a radiation heat- 
transfer model. A six-flux model of radiative heat transfer Is described. 

This model provides only a limited geometric description of the radiation 
transfer process, but a comparison with experimental data Indicates an encour- 
aging level of agreement. 


NUMERICAL ACCURACY 

The following Is a general form of the time-averaged equations that must 
be solved In a combustor flow code: 

f^(pU<p) + f^(pV<p) - fj) - fj) = (1) 

where <p can represent U, V, W, u^Uj, K, e, or H; where Is the 

diffusion coefficients (e.g., y e ff); and s <p the source term (e.g., - 3P/3x). 
The equation represents convection of a conserved scalar that Is subtracted by 
diffusion terms equal to a source term that could represent a pressure gradi- 
ent or a source/sink term. These equations are solved by discretizing on a 
staggered mesh using the finite volume method (fig. 1). The staggered mesh 
system Is used to avoid the pressure-velocity decoupling that can result In a 
finite volume representation of Incompressible flow. 

An example of the discretization of a convective term using central dif- 
ferencing Is the following: 
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(Alternate types of differencing are possible and some will be Illustrated In 
the following paragraphs.) Once all the terms have been approximated, like 
terms are then rearranged In a substitution formula which can be solved using 
a tri-diagonal matrix algorithm (TOMA) procedure applied In alternate 
directions (ADI): 

*p<Pp « a N «p N + a s «p s ♦ a E q> E + <ypg ♦ S (3) 

where, for example. 


a N = 




The source and diffusive terms In these equations are approximated using 
central differencing which Is second-order accurate (Error « 0(ax)2). The 
convective terms have typically been differenced using hybrid differencing 
which reverts to upwind differencing when the absolute value of the cell 
Reynolds number, or cell Peclet number. Is greater than 2, which Is the main 
reason for the loss of numerical accuracy. An example of convective discreti- 
zation using hybrid differencing Is as follows: 


au<p 

dx 




(U<p) - ( u<p) 

e w 
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For upwind, with u > 0, 


( u«p) = u <p_ 
e e P 


and 


<“*>„ - 2 u w<*w * V 

Since upwind differencing Is only first-order accurate. It can Introduce an 
extensive amount of numerical diffusion Into the calculation. The virtue of 
this scheme Is that It provides "bounded," stable solutions. Higher order 
convective schemes Invariably yield nonphysical oscillations In the solution. 
Upwind differencing avoids this, at a considerable loss In accuracy. 

To alleviate this problem, NASA has conducted a program to Identify and 
Incorporate an Improved accuracy differencing scheme Into a combustor flow 
code. Under a portion of this program a variety of differencing schemes were 
examined In several test calculations. The schemes examined Included QUICK 
(quadratic upstream Interpolation) and SUD (skewed upwind differencing). 

QUICK differencing was developed by Leonard (ref. 1). This scheme Improves 
the accuracy of convective differencing by performing an upwind biased quad- 
ratic Interpolation. For u > 0, 
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(U») e . u e (- 3 » w * f„p . §* E ) 

<“*>« - u w (- e^ww * «*w * 8»p) 

where grid point locations are as noted on figure 2. This scheme Is second- 
order accurate and can produce nonphysical oscillations In the solution. SUO 
(skewed upwind differencing) (ref. 2) attains high accuracy by differencing In 
an upwind manner along the flow streamlines. While maintaining the same for- 
mal accuracy as upwind differencing, the truncation error In SUD Is smaller. 

For example: 

For u > 0 and V > 0, 

(u<p) e = VP | 

(u <p) w = u w O " a)<p w + u w a<p sw ( (6) 

V = minimum of (l,V/2u) I 

where grid point locations are as noted In figure 3. As with QUICK, SUD can 
produce nonphysical oscillations In the solution; therefore, a scheme to 
"bound" SUD was also examined. This scheme employs the concept of flux- 
blending (ref. 3), wherein a bounded flux determined from upwind differencing 
Is blended with the unbounded, but more accurate, SUD flux. The main factor 
Is to blend as little of the lesser accurate scheme while still maintaining a 
properly "bounded" solution. This procedure, called BSUDS, starts from an 
Initial, totally skew-differenced estimate and blends an upwind flux If the 
solution Is out of the range of neighboring values. If the solution Is In 
range (l.e., bounded), then no blending Is performed. 

An Illustration of the accuracy of the upwind, QUICK, and SUD schemes Is 
seen In figure 4. This figure displays the results of a single-point scalar- 
transport calculation made for various flow angles. All schemes agree with an 
exact solution (no error) at a zero flow angle; however, at angles greater 
than zero, each scheme displays some degree of error relating to numerical 
diffusion. The error displayed by upwind differencing Increases with flow 
angle to a maximum at 45°. The QUICK scheme displays a similar behavior, but 
the overall error Is much less. The SUD scheme displays a maximum error 
around 15°, but It tends to zero at angles approaching 45°. Both QUICK and 
SUD display a much higher level of accuracy than upwind. 

Although a scalar transport calculation Is useful for a general examina- 
tion of some aspects of differencing scheme performance, a laminar flow calcu- 
lation Is a more complete test. The results of a series of laminar flow 
calculations from reference 4 are displayed In figure 5. In this figure, 
axial velocity profiles at a distance of one-half a duct height from the Inlet 
are shown for two different computational meshes. In these calculations, the 
steepness of the velocity profile Indicates accuracy. Steep velocity profiles 
are exhibited by QUICK and BSUDS; the upwind profiles exhibit a high degree of 
numerical diffusion. On the coarse mesh, BSUDS appears to be more accurate 
than QUICK; on the fine mesh, there Is not much of a distinction between the 
two schemes. 
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An Important aspect of Improved accuracy Is the effect that the differ- 
encing scheme has on the rate of convergence. The computational times to 
converge the system of governing equations In the previous laminar calcula- 
tions are shown In table I. The convergence times are ratloed to the upwind 
convergence times to clearly Illustrate the computational penalty paid to 
attain Improved accuracy. Generally, the Improved accuracy schemes required 
from 3 to 15 times longer to reach a converged solution. To a degree one can 
hope that this computational penalty can be offset by using coarse meshes to 
achieve the same overall level of accuracy. In reality, a relatively fine 
mesh Is needed, even with the high-accuracy schemes. In any case, the need 
for Improved solution algorithms for these more accurate differencing schemes 
Is strongly Indicated. 


SOLUTION ALGORITHMS 

The previous section demonstrated the need for Improved solution algo- 
rithms. This section examines the widely used SIMPLE algorithm and two modi- 
fications to this scheme. These schemes cover only a small portion of the 
wide range of methods to accelerate solution convergence. Vectorlzatlon, 
direct-solution methods, and multigrid methods - to list just a few - are all 
areas of active research that are certain to yield much greater computational 
benefits In the near future. 

In the SIMPLE algorithm, a guessed pressure field Is Inserted Into the 
discretized momentum equations to obtain a velocity field. The pressure field 
Is corrected by an equation which Is derived through a combination of the con- 
tinuity and momentum equations. The velocity field Is then updated and Is 
used In the solution of the equations for k, e, and <p. The corrected pres- 
sure field Is treated as the guessed pressure field, and the procedure Is 
repeated until a converged solution Is obtained. 

The following velocity correction equation Is used In the SIMPLE scheme 
for u at point e: 

A u' =5" A . U' + ( P * - P ' ) A 
e e ^ nb nb P E e (7, 

Primes Indicate corrections to old values. The underlined term, which repre- 
sents the Influence of corrected pressures on neighboring velocities. Is 
neglected In the SIMPLE algorithm. The converged solution Is unchanged by the 
exclusion of this term since, for the steady-state solution, the corrections 
go to zero. Neglecting terms, however, does force the use of low underrelaxa- 
tlon factors, which can slow convergence. 

The SIMPLER algorithm Improves on the SIMPLE scheme by Including the 
previously neglected terms when calculating the pressure field. The calcula- 
tion sequence starts with a guessed velocity field. An equation that solves 
for the pressure field (using the terms Ignored In the SIMPLE scheme) Is cal- 
culated from the guessed velocity field. This pressure field Is then used to 
solve the discretized momentum equations to obtain a velocity field. The 
velocity field Is corrected In a manner similar to the SIMPLE velocity correc- 
tion. This velocity field Is then treated as the guessed velocity field, and 
the Iteration procedure Is repeated until convergence Is reached. 
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Because additional equations are solved, each Iteration through the 
SIMPLER routine Involves more computational time than an Iteration step through 
SIMPLE. However, higher underrelaxatlon factors can be applied In the SIMPLER 
routine, thereby accelerating convergence. 

The PISO scheme also takes Into account the terms neglected In the SIMPLE 
code, but In a different manner. The PISO routine mimics the SIMPLE approach 
until the end of the first Iteration. At this point, the PISO scheme employs 
an equation containing the neglected terms to correct the pressure and veloc- 
ity field to more closely agree with continuity. Again, this procedure Is 
repeated until the solution converges. In this manner, the PISO code allows 
for higher underrelaxatlon factors to accelerate convergence. 

The performance of these various solution schemes Is displayed In 
figure 6. As an example, the convergence times for a 38-by-38 grid point 
calculation are plotted as a function of underrelaxatlon factor (fig. 6(a)). 
When a reasonably large underrelaxatlon factor Is used, SIMPLER or PISO con- 
verge about twice as fast as SIMPLE. In addition, SIMPLER and PISO converge 
over a larger range of underrelaxatlon factors than SIMPLE. From an engineer- 
ing standpoint, Improving the "robustness" of a computational scheme Is often 
just as Important as accelerating convergence. The computational benefit of 
using SIMPLER and PISO for fine-mesh calculations Is displayed In figure 6(b). 
The greater the number of mesh points used In the calculation, the greater the 
benefit of SIMPLER or PISO over SIMPLE. For example, a calculation of approx- 
imately 3300 grid points converges three times faster using PISO or SIMPLER 
than It does using SIMPLE. This Is a savings of about 600 CPU seconds. The 
1440 grid points calculation demonstrated a savings of only about 100 CPU 
seconds . 


INLET BOUNDARY CONDITIONS 

In any calculation of a complex, three-dimensional turbulent flow, the 
boundary conditions that are needed In the calculation are frequently unknown 
or unmeasured. The use of Inappropriate values at the computational boundary 
can sometimes be the main limit to the calculation's predictive capability. 
This error can sometimes be more Important than numerical accuracy or turbu- 
lence model considerations. Figure 7 shows an example of a three-dimensional 
jet-ln-crossf low calculation using alternate boundary conditions. In one 
calculation, the jet orifice flow was specified as having a uniform plug flow 
at a position two jet diameters upstream of the orifice outlet. This allowed 
the flow to distort as It exited the orifice outlet. The second calculation 
specified a uniform plug flow at the orifice outlet. The resulting axial 
velocity profiles are compared with experimental data In figure 8. Somewhat 
surprisingly, the uniform boundary condition at the orifice compared more 
favorably with experimental data than the theoretically more correct distorted 
profile. This should not be Interpreted as an endorsement of the use of 
uniform-plug-flow boundary conditions for these types of flows. Indeed, there 
Is some Indication of experimental error (ref. 5). This example Is meant only 
to Illustrate that unmeasured or unknown boundary conditions can significantly 
affect a numerical calculation. 


TURBULENCE MODELS 

Although a great deal of progress Is being made In solving the three- 
dimensional, time-dependent Navler-Stokes equations In large-eddy or direct 
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numerical simulations, practical engineering calculations currently require 
the introduction of some form of turbulence modeling. These models are based 
on either Reynolds or Favre averaging of the exact Navier-Stokes equations, 
reducing the unsteady form of these equations to an averaged form. Currently, 
the most widely used turbulence model is the two equation, k-e, closure. This 
model relates the Reynolds stresses to a turbulent viscosity through 
Boussinesq's eddy-viscosity concept: 


-pu 1 u j = pv 


faUj 
T\ixj 



( 8 ) 


where vj = C u k 2 /e. The turbulent viscosity is related to the kinetic 

energy k and the dissipation rate c of the turbulence. The transport 
equations are then solved for k and e: 
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Where C = 0.09, C _ = 1.44, C _ = 1.92, a. = 1.0, o = 1.3, a = 0.9, k 

v el cl K e <p 
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(9) 


The model constants typically employed are those recommended in 
reference 6. The turbulent Schmidt number a ip , is frequently changed from 
0.9 to as low as 0.2, depending on the flow being studied. 


This two-equation model is based on several assumptions which should be 
considered when making a numerical calculation. First, the flow is assumed to 
be close to equilibrium; that is, the flow properties change relatively slowly. 
Second, the turbulence Reynolds number is assumed to be high. Third, the tur- 
bulence is assumed to be isotropic. 


The main concern is how well this model, with its Inherent assumptions, 
can represent combustorlike flow fields. Figure 9 displays a comparison 
between laboratory experiments and numerical predictions of two different iso- 
thermal flows. (Figs. 10 and 11 show the locations of the measurements that 
were made in the flow fields.) In the two-dimensional bluff-body comparison, 
a major disagreement between measurements and predictions is evident at an 
axial distance x/D of approximately 0.8. This corresponds to the end of the 
recirculation zone and causes an Incorrect prediction of recirculation zone 
length. The jet-in-crossf low comparison displays a similar disagreement in 
the region where the turbulence intensity is high. In this comparison, both 
hybrid and BSUDS differencing were used in the predictions. It is obvious 
from the displayed results that numerical accuracy can have a major Impact on 
the comparison with experimental measurements. Hybrid differencing is so 
completely Influenced by numerical diffusion that the qualitative agreement 
between experiment and calculation, evident in the BSUDS results, is elimi- 
nated. The BSUDS results are not grid-independent, but it seems unlikely that 
this will fully explain the noted disparity. 
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Both of these flow fields display the greatest disagreement between exper- 
iment and calculation where the turbulence Intensities are the highest. These 
are regions where the flow field Is likely to be far from equilibrium. It Is 
Interesting to note that a full Reynolds stress transport (RST) model calcula- 
tion (presented by McGuIrk, J.3., Papadlmltrlou, C., and Taylor, A.M.K.P. at 
the Fifth Symposium on Turbulent Shear Flows held at Cornell University, 

Ithaca, New York, August, 1985) did not yield appreciably better results for a 
similar flow field. Because both models appear to lose validity around the 
region of the stagnation point, further model development Is needed. 

Reaction closures Involve a further series of assumptions and approxima- 
tions. The simplest level of closure Is to assume that the reaction Is mixing 
limited and Ignore the effect of temporal density fluctuations. The reaction 
rate can then be related to an eddy lifetime (e/k) using either a Magnussen- 
Hjertager or Spaldings eddy-breakup model: 

Rate = -S (10) 

m p 

where 

T Mo 2 M r 1 

S m p = " p k m1n | Am F» A STOIC ’ AB STOIC + 1J 


or 


S m p “ " C R p k m Fu 

where A, B, and Cp are empirical constants and e/k can be considered as 
the eddy lifetime. The density and temperature throughout the flow field can 
then be established from equilibrium chemistry. 

Improved physical realism can be added to this model by Introducing the 
probability density function (PDF) for mixture fraction f. This can account 
for the unsteady time history of the mixture fraction at each point In space. 

A nonlinear functional dependence between concentration, temperature, and 
density occurs as pockets of alternately hot and cold gases are swept past 
each point In space. Integrating the resultant PDFs allows a determination of 
the mean properties of the flow. 

Figure 12 displays the results of both an eddy-breakup and a PDF model 
calculation of a planar-mixing layer. The PDF model calculations (ref. 7) 
demonstrate a significant Improvement over the eddy-breakup results. The main 
factor Involved In this Improvement Is the more physically correct representa- 
tion of the mixture fraction variation at the point of maximum temperature 
rise. Of course, both of these models are only as good as the turbulence 
closure and are wholly Inadequate when finite-rate chemistry Is Important. 

RADIATION HEAT TRANSFER 

The final topic to be covered In this review concerns radiative heat 
transfer. Radiation Is one of the most significant and least understood heat 
loads to the combustor liner. It can also play a significant role In the 
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determination of flame temperature. Any numerical description of a gas tur- 
bine combustor must Include a radiation heat-transfer model. The six-flux 
model Is the model most commonly used to approximate multidimensional radi- 
ative transfer. In this model, differential equations describing the radi- 
ative fluxes In positive and negative directions along the principal axis are 
solved: 
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where 

RX,r,z composite fluxes 
a absorption coefficient 

S scattering coefficient 

E oT 4 

The main Input to this analysis concerns the optical characteristics of the 
hot gas and soot which must be arbitrarily specified or calculated through a 
soot formation and oxidation model. 

The performance of the six-flux model (ref. 8) Is displayed In figure 13. 
Although the model overestimates the radiative heat transfer In comparison 
with experimental data, the qualitative trend Is quite closely followed. Given 
the large number of approximations used In the analysis, the agreement with 
experimental data Is quite surprising. The six-flux model does not accurately 
treat the angular dependence of energy transfer, and the determination of the 
optical characteristics of the soot cloud still remains as an area of needed 
research; however, fairly good results appear to be possible In this example. 

CONCLUDING REMARKS 

Perhaps the most Important question any review on numerical modeling can 
answer Is whether or not current computational codes can be usefully employed 
In the design of combustion devices. Certainly a great deal of research Is 
needed before one can expect quantitative predictive accuracy, and It seems 
likely that some hardware problems will only be resolved through development 
testing. The best computer program will never replace the designer's Innova- 
tive mind, but computer predictions can be used to extend the designer's pro- 
ductivity. New designs can be examined much more rapidly on the computer than 
In hardware testing. Development costs can be reduced. The promise of this 
computer-based design methodology Is so great that these numerical models will 
be used despite their deficiencies. Designers should not and probably will 
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not abandon empirical design tools, but the cautious adoption of numerical 
models In the design process Is a trend which can reap Important benefits. 
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TABLE I. - RATIO OF CONVERGENCE TIMES 
FOR VARIOUS DIFFERENCING SCHEMES 
WITH UPWIND CONVERGENCE TIMES 
USED AS THE STANDARD 


Mesh 

Upwind 

BSUDS 

QUICK 

Coarse (30 by 22) 

1 

6.4 

A 

0 . c 

Fine (58 by 38) 

1 

14.7 

15.7 
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Figure 1. - Staggered mesh system for discretizing equations 
by finite volume method. 
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Figure 2. - Staggered mesh system for 
quadratic upstream interpolation. 



Figure 3. - Staggered mesh sys- 
tem for bounded skew upwind 
differencing scheme (BSUDS). 
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Figure 6. - Convergence times. 
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(b) For various computational meshes. 
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Figure 7. - Three-dimensional jet-in-crossflow 
calculation. 
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Figure 8, - Calculated axial velocity profiles com- 
pared with experimental data. 
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(a) Two-dimensional bluff-body wake. (b) Three-dimensional jets in crossflow. 

Figure 9. - Comparison of laboratory experiments with numerical predictions of two different isothermal flows. 



0 1.0 
X/D 


Figure 10. - Two-dimensional bluff-body geometry. 
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(b) Eddy breakup model. 


(c) PDF model. 


Figure 12. - Results of calculations of a planar-mixing layer. 
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(a) Schematic drawing of experimental measurement 
locations. 



(b) Performance of six-flux 
model. 

Figure 13. - Radiative heat transfer test calculation. 
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